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We construct a coarse-grained (CG) model for dipalmitoylphosphatidylchioline (DPPC) / cholesterol bilayers 
and apply it to large-scale simulation studies of lipid membranes. Our CG model is a two-dimensional repre- 
sentation of the membrane, where the individual lipid and sterol molecules are described by point-like particles. 
The effective intermolecular interactions used in the model are systematically derived from detailed atomic-scale 
molecular dynamics simulations using the Inverse Monte Carlo technique, which guarantees that the radial dis- 
tribution properties of the CG model are consistent with those given by the corresponding atomistic system. 
We find that the coarse-grained model for the DPPC / cholesterol bilayer is substantially more efficient than 
atomistic models, providing a speed-up of approximately eight orders of magnitude. The results are in favor of 
formation of cholesterol-rich and cholesterol-poor domains at intermediate cholesterol concentrations, in agree- 
ment with the experimental phase diagram of the system. We also explore the limits of the novel coarse-grained 
model, and discuss the general validity and applicability of the present approach. 



I. INTRODUCTION 

Cell membranes are remarkably flexible and durable struc- 
tures enclosing and protecting the contents of cells or 
organelles4i2i2iii5i& They are, however, by no means mere cas- 
ings, but bustling hubs, where signaling, recognition, and 
transport take place. The fact that a huge variety of cellular 
processes are governed by membranes makes them a fasci- 
nating and biologically relevant example of soft and distinctly 
thin interfaces. 

The complexity and biological relevance of membranes is 
largely due to the variety of proteins and lipids that are their 
main building blocks. It is intriguing that there are typically 
more than a hundred different lipid species in a given type 
of biological membrane, all assumed to have some particular 
purpose.'-'' Instead of being static, membranes are highly dy- 
namic and characterized by distinct phases and internal fluctu- 
ating structures, thus allowing proteins to function under non- 
equilibrium conditions.^ Understanding the overall properties 
of membranes is therefore a major challenge involving studies 
over large scales in time and space: starting from the atom- 
istic and molecular regimes where small-scale processes such 
as ion flow through channel proteins takes place all the way to 
mesoscopic and macroscopic regimes where large-scale pro- 
cesses such as phase separation and membrane fusion are im- 
portant. 

The present understanding of membrane systems is largely 
based on experimental studies, where techniques such as flu- 
orescence spectroscopy, nuclear magnetic resonance, and var- 
ious scattering methods have been used.''^'^"^'^ At the same 
time, experiments have been complemented by theoretical and 
computational modeling ]^i-^-^i^i^i'"i"i'^i'-^ Thanks to the inter- 



play between the two fields, a more detailed understanding of 
membranes and their biological relevance is emerging. 

As for computational modeling of membrane systems in the 
atomistic regime, classical molecular dynamics (MD) is the 
method of choice i2i2ii2iiiii2ii^ It provides detailed information 
on the structure and dynamics of individual lipid molecules, 
as well as insight into processes such as the complexation and 
hydrogen bonding properties of different lipid species or the 
effect of enzymes on membranes. The main limitation of the 
atomistic approach is the computational load. With currently 
available computer power, the standard size for systems is 128 
lipid molecules, coiTesponding to a linear system size of about 
5-7 nm in the bilayer plane. The duration of such a simula- 
tion is typically limited to about 100 ns. A few more ambitious 
MD simulations on systems with a larger number of molecules 
have been reportedp^ii^ii^ but the sizes are still rather modest: 
the largest MD studies of lipid bilayers contain of the order of 
lO'^ lipid molecules. The time scales reached in such simula- 
tions are cuiTently only tens of nanoseconds. 

The above limitations are problematic since many inter- 
esting phenomena in lipid membrane systems occur at much 
longer time and length scales. Examples of such phenomena 
are domain formation, bilayer fusion, and cooperative mo- 
tions associated with phase changes. Domain formation is 
a particularly interesting issue, since there is plenty of ex- 
perimental evidence pointing to the formation of lateral do- 
mains in many-component bilayers. ^'^ The most topical is- 
sue are lipid rafts fi^ti^ii^^ which are thought to be dynamic, 
ordered lateral domains comprised mainly of phosphatidyl- 
choline, cholesterol, and sphingomyelin molecules. Rafts 
have been suggested to be involved in a wide range of cellu- 
lar processes including membrane trafficking and sorting of 



proteins, '^<ii^ii2i2£ which emphasizes the need to understand 
large-scale properties of membrane domains. The dimensions 
of these dynamic domains are believed to range from tens to 
hundreds of nanometers, still beyond the limits of atomistic 
simulations. To reach the necessary length scales, i.e., hun- 
dreds of nanometers, we therefore need to resort to coarse- 
grained models that employ effective interaction potentials 
for simplified molecular descriptions. ^•^'■^^•^^ To gain insight 
into large-scale properties of lipid membrane systems, the 
main objective is hence to develop and employ coarse-grained 
membrane models incorporating only the essential properties 
of the underlying system. 

Previous studies in this field have been few, although re- 
cent progress is very promising. As for semi-atomistic mod- 
els of bilayers, a number of interesting approaches have 
been suggested , ^"^•^^■^''■^^•^^•^'^-•^"•'^'•'^^••^'' The guiding principle 
is that small groups of atoms are represented as single in- 
teraction sites, thus reducing the computational complexity 
of the model. In the model by Marrink et al.^'^ there are 
Lennard-Jones, harmonic, and electrostatic interactions be- 
tween the coarse-grained particles, and the interaction param- 
eters have been tuned to match experimental quantities such 
as heats of vaporization or densities. The systems may be 
simulated using classical MD. The approach by Lipowsky 
et al. is somewhat more phenomenological but similar in 
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Shelley et al, in turn, have employed a large 
variety of different interactions between the coarse-grained 
particlespii^2iS and the interaction parameters have been ad- 
justed to match results from both experiments and atomic sim- 
ulations. Groot and Rabone have employed the dissipative 
particle dynamics (DPD) technique and chosen soft potentials 
between the coarse-grained particles.^** The repulsive interac- 
tion parameters have been derived from compressibility and 
solubilities. The DPD studies by Groot and Rabone have been 
complemented by the simulations of Smit et al^^^ Ayton et 
al. have also employed the DPD technique. They have used 
material properties from atomistic simulations to parameter- 
ize meso-scale and macro-scale models of lipid bilayers and 
unilamellar vesiclesi^ii^^i^^iS 

Alternatively, one can design phenomenological models 
where the number of degrees of freedom is as small as pos- 
sible. One of the best examples of this approach is the work 



by Mouritsen et al. 
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They have developed and used 



an off-lattice model where lipid and sterol molecules are de- 
scribed as hard-core particles with internal (spin-type) degrees 
of freedom. This approach has allowed them to design mod- 
els whose phase diagrams are in qualitative agreement with 
experimental ones for phosphatidylcholine (PC) /cholesterol 



and PC / lanosterol systems 
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Additionally, the mod- 



els have been successful in describing lateral diffusion in 
PC /sterol bilayer mixtures. ^^ The work by Mouritsen et al. 
demonstrates that purely phenomenological models can be 
very useful in accessing scales larger than those within reach 
of atomistic simulation techniques. On the other hand, due 
to their phenomenological nature, the scope of these models 
may be limited. 

In general, there is no single method of constructing models 
for mesoscopic or macroscopic phenomena, and each case has 



to be considered separately. Hence, systematic and general 
approaches that simplify the construction of coarse-grained 
models and reduce the number of phenomenological and tun- 
able parameters would be of great interest. A generally useful 
approach be easily extendable or modifiable to describe dif- 
ferent kinds of systems. 

A promising candidate is the Inverse Monte Carlo tech- 
nique (IMC)."^^"*^ It allows one to derive all effective interac- 
tion potentials systematically from atomic-level information 
such that the most relevant structural properties of the atomic- 
level system are reproduced by the coarse-grained model. 
This approach has been used in other soft matter systems. For 
example Lyubartsev et al.'^^-'^'^-'^^ have used IMC for construct- 
ing a coarse-grained model for sodium and chloride ions in 
water, and for describing the binding of different alkali ions 
to DNA.^^ The approach employed by Shelley et al. is also, in 
part, based on the central ideas of the IMC method.^' Notably, 
IMC is not only systematic, but also highly adjustable, as the 
level of coarse-graining and the number of degrees of freedom 
can be tuned. It is thus possible to choose how large scales are 
to be studied as well as how much detail is to be included. 

In the present study, we apply the IMC approach'*^'*'* to con- 
struct a coarse-grained (CG) model for a lipid bilayer contain- 
ing dipalmitoylphosphatidylcholine (DPPC) and cholesterol. 
This system was chosen because DPPC is one of the most 
studied phospholipids, and cholesterol is the most important 
sterol molecule found in plasma membranes of eukaryotic 
cells. Further, the system has a rich and interesting phase 
behavior^^ characterized by three main phases (see Fig. [Q. 
It has been suggested that at certain temperatures and choles- 
terol concentrations the system might display cholesterol-rich 
domains'*** or superlattice domains.'*^ 

We first performed extensive atomic-level MD simulations 
of the DPPC /cholesterol system at six different cholesterol 
concentrations.^" The results of these simulations agree well 
with experiments and other simulations. Based on these atom- 
istic considerations, we construct a coarse-grained model in 
which we describe each molecule by a single point-like parti- 
cle moving in a two-dimensional plane. We use the IMC tech- 
nique to derive effective interaction potentials for the coarse- 
grained particles. These interactions are constructed such 
that the CG model reproduces the radial distribution functions 
(RDFs) calculated from the atomic-level MD simulations. Be- 
cause the RDFs can be used to characterize the phase behav- 
ior, the model should also qualitatively reproduce the phase 
behavior of the microscopic model. 

Using the coarse-grained model, we study 
DPPC /cholesterol bilayers with cholesterol concentra- 
tions varying from 0% to 50%. We first validate the CG 
model by comparing its behavior to that of the atomic-scale 
model. As the degree of coarse-graining is very high, the 
model allows us to study the properties of the bilayer on 
length scales of the order of lOOnm along the plane of the 
membrane. The computational gain can be approximated 
to be around eight orders of magnitude compared to the 
atomistic MD case. 

We find that the coarse-graining approach can pro- 
vide plenty of insight into large-scale properties of many- 
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FIG. 1: Sketch of experimental phase diagram for the 

DPPC / cholesterol system.'*'' At high temperatures and low choles- 
terol concentrations, there is a liquid-disordered (Id) phase, which is 
a fluid-phase characterized by lipid acyl chains with a high degree 
of conformational disorder. When the temperature is lowered, the 
system goes through the main phase transition at Tm ~ 311 K to a 
solid-ordered (so) phase. The so phase is essentially a solid phase in 
which acyl chains are conformationally ordered and the positions of 
the molecules are characterized by translational order in the bilayer 
plane. Finally, at high cholesterol concentrations, there is a liquid- 
ordered (lo) phase, characterized both by a high degree of acyl chain 
ordering and the lack of translational order found in the Id phase. 
At intermediate cholesterol concentrations there are wide Id-lo and 
so-lo coexistence regions. The dots represent the concentrations at 
which the atomic-scale molecular dynamics simulations have been 
performed. An additional MD simulation was performed at 50% 
cholesterol. 



component membrane systems. In this case it allows us to ob- 
serve formation of cholesterol-rich and cholesterol-poor do- 
mains at intermediate cholesterol concentrations, in agree- 
ment with the experimental phase diagram of the system4i 
We also explore the limitations of the model, and discuss its 
general validity as well as possible future applications. 



II. MOLECULAR DYNAMICS SIMULATION DETAILS 

The underlying MD simulations have been described in 
detail elsewhere,^" and only a brief summary is given here. 
We simulated fully hydrated lipid bilayer systems consisting 
of 128 macromolecules, i. e., DPPCs and cholesterols, and 
3655 water molecules. The simulations were performed at 
six cholesterol molar fractions: 0%, 4.7 %, 12.5 %, 20.3 %, 
29.7%, and 50%. These concentrations are indicated in the 
phase diagram of Fig. [J The duration of each simulation was 
100 ns and the linear sizes of the systems in the plane of the 
bilayer were between 5 and 7 nm. 

The starting point for the simulations was a united atom 
model for a fully hydrated pure DPPC bilayer that has been 
validated previouslyi^ii^SiS The parameters for bonded and 
non-bonded interactions for DPPC molecules were taken from 



a study of a pure DPPC bilayer,^"* and partial charges from the 
underlying model description.^^ The cholesterol force field 
was taken from an earlier study.^^ 

As an initial configuration for the pure DPPC bilayer we 
used the final structure of run E discussed in Ref. ,53, For 
systems containing cholesterol, the initial configurations were 
constructed by replacing randomly selected DPPC molecules 
with cholesterols. The same number of DPPC molecules was 
replaced in each of the two monolayers. To fill the small 
voids left by replacing DPPC molecules by somewhat smaller 
cholesterol molecules, the system was equilibrated in several 
stages.^" 

The MD simulations were performed at a temperature T — 
323 K using the GROMACS molecular simulation packagei^ 
The main phase transition temperature for a pure DPPC bi- 
layer is Tm ~ 311 K,^^ indicating that the MD simula- 
tions have been conducted above Tm. The time step for 
the simulations was chosen to be 2.0 fs. Long-range elec- 
trostatic interactions were handled using the Particle-Mesh 
Ewald method. ^^ After the initial equilibration we performed 
100 ns of MD in the NpT ensemble with a Berendsen thermo- 
stat and barostat^^ for each cholesterol concentration. For all 
systems up to and including the cholesterol concentration of 
29.7 %, a simulation lasting 100 ns guarantees good sampling 
of the phase space. The results for 50 % cholesterol should be 
regarded with some caution as the diffusion of the DPPC and 
cholesterol molecules is already quite sloW(ifi 



III. COARSE-GRAINED MODEL 

Using the MD simulations as a basis, we have constructed a 
coarse-grained model for a DPPC / cholesterol bilayer Since 
the main goal of the present project is to study the large- 
scale structural properties of the bilayer, the degree of coarse- 
graining must be high. A way to achieve this goal is to de- 
scribe each DPPC and cholesterol molecule by its center-of- 
mass (CM) position. The macromolecules are taken to be sin- 
gle point-like particles that move and interact in two dimen- 
sions with continuous coordinates. At the same time, the sol- 
vent degrees of freedom have been integrated out altogether, 
i. e., the model contains no explicit water. 

Let us briefly list the assumptions we have made in con- 
structing the CG model. To start with, we consider a lipid 
bilayer as a purely two-dimensional sheet comprised of two 
weakly interacting leaflets. For this reason, we focus on one 
monolayer only. This assumption is well justified since in- 
terdigitation in DPPC / cholesterol bilayers is minor^*^ Conse- 
quently, the friction between the two leaflets is weak and they 
can be regarded as largely independent from each other Fur- 
thermore, we neglect the out-of -plane fluctuations of the bi- 
layer and assume it to be strictly planar Such fluctuations de- 
crease when the cholesterol concentration is increased,'^ mak- 
ing this a reasonable assumption especially at higher choles- 
terol concentrations. 

Due to its coarse grained nature, our model is dissipative. 
This stems mainly from the fact that the water molecules are 
not included in the CG model. Further, in constructing the 



model the conformational degrees of freedom of the macro- 
molecules have been integrated out. Yet another reason is that 
we consider one of the leaflets rather than the whole mem- 
brane. If we were to study dynamical phenomena, the dy- 
namics should be chosen such that there are both stochastic 
and dissipative force components describing those degrees of 
freedom that have been excluded from the CG model. As we 
will focus on structural quantities of the membrane system, 
we will not have to worry about the choice of dynamics. In- 
stead, we use the Metropolis Monte Carlo (MC) techniqueiS 
The question of implementing realistic dynamics to the model 
is considered in more detail at the end of Sect. lVIl 

As for the interactions between the point-like DPPC and 
cholesterol particles, we assume they can be adequately de- 
scribed using pairwise, radially symmetric effective poten- 
tials. The effective interactions are computed as follows. 
From the atomistic MD simulations we calculate radial dis- 
tribution functions for the CM positions of the molecules. To 
link our coarse-grained model to the atomic-level system, we 
require that the coarse-grained model accurately reproduces 
these RDFs. This is accomplished by constructing the ef- 
fective interaction potentials using the Inverse Monte Carlo 
method."'^"*'* In principle, also other canonical averages than 
the RDFs could be used as an input. However, the RDFs 
calculated from the atomic-scale MD simulations are easy to 
compute, and more importantly, give a detailed structural de- 
scription of the system in the plane of the bilayer at short 
length scales. Because the RDFs can be used for character- 
izing the phase behavior of the system, the coarse-grained 
model should at least qualitatively reproduce the phase be- 
havior of the original atomic-scale system. 

In addition to the above, we fix the area per molecule in the 
CG model to be the same as the average area per molecule 
calculated from the atomistic MD simulations. The MC simu- 
lations will therefore be conducted in the canonical ensemble. 



IV. MODEL CONSTRUCTION AND VALIDATION 

A. Radial Distribution Functions and Areas per Molecule 
from Atomistic MD Simulations 

For the construction of the coarse-grained model, we need 
to obtain the radial distribution functions for the center-of- 
mass positions of the molecules. Figure 121 shows the RDFs 
calculated from the atomic-scale MD simulations. Before 
calculating the RDFs, the CM positions have been projected 
to the plane of the bilayer The two monolayers have been 
considered separately and the resulting RDFs are averages 
of the two. The first 20 ns of the MD data have been 
discarded to allow the system to equilibrate fullyi^ After 
20 ns the area per molecule has converged for all cholesterol 
concentrations,^" and the radial distribution functions show 
no systematic changes. The radial distribution functions were 
calculated up to 2.5 nm for concentrations lower than 50 %. At 
the highest concentration of 50 %, the linear size of the system 
in the bilayer plane is occasionally below 5 nm, and therefore 
in this case the RDFs were cut off at 2.4 nm. The errors of 
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FIG. 2: Radial distribution functions calculated from MD simula- 
tions for (a) DPPC-DPPC, (b) DPPC-cholesterol and (c) cholesterol- 
cholesterol pairs. The RDFs are calculated from the center-of-mass 
positions of the molecules, which have been projected to the plane of 
the bilayer. 



the RDFs can be estimated to be of the order of a few percent, 
with somewhat higher errors at low cholesterol concentrations 
for the RDFs involving cholesterol. To minimize the effect of 
random errors on the Inverse Monte Carlo procedure, we used 
a spline-fitting procedure designed for noisy data^' to smooth 
the RDFs. 

For all concentrations the RDFs indicate liquid-like behav- 
ior. At short length scales there are broad peaks and at larger 
r the functions approach unity. In other words, although there 
is short-range order, there are no signs of of long-range or- 
der characteristic to solid-like phases. This confirms that we 
are probing the region of the phase diagram where the sys- 
tem is in the Id, lo, or coexisting Id and lo phases (see the 
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FIG. 3: Average area per molecule as a function of cholesterol con- 
centration, calculated from MD simulationsi— 



points marked in Fig. [Q. As the cholesterol concentration 
is increased, the radial distribution functions change more or 
less systematically. The peaks in the DPPC-DPPC distribution 
become sharper, manifesting an increase in the lateral short- 
range ordering. Further, more peaks appear at larger r, which 
means that the range of the ordering increases slightly. Similar 
effects are observed in the case of the cholesterol-cholesterol 
RDFs, although the 4.7 % concentration deviates somewhat 
from the general trend. For the DPPC-cholesterol distribution 
the changes are not quite as systematic. Nonetheless, when the 
cholesterol concentration grows, the range of ordering seems 
to increase slightly and the peaks generally become sharper. 

A notable feature of the RDFs, especially for the DPPC- 
cholesterol pairs, is the fact that some RDFs do not approach 
zero at the origin. This is because the RDFs have been cal- 
culated for the CM positions of the molecules, which have 
been projected to the plane of the bilayer. It is not too difficult 
to imagine a situation where the projected CM positions of a 
rigid, short cholesterol molecule and a DPPC molecule with 
long, flexible tails are essentially on top of each other. 

In addition to the RDFs, we need the average area per 
molecule from the MD simulations to fix the area per 
molecule in the coarse-grained model. The average areas 
per molecule at different cholesterol concentrations have been 
published previously, ^° and are shown in Fig.|3]for reference. 
The area per molecule for a given configuration has been com- 
puted by dividing the size of the simulation box in the bi- 
layer plane by the number of molecules in one monolayer. 
The area per molecule decreases monotonically with an in- 
creasing cholesterol concentration, in agreement with other 
simulation g'^i^^ and related experimentSiS 



B. Constructing Effective Potentials 

Based on the RDFs, we have constructed effective interac- 
tion potentials for our coarse-grained particles using the IMC 
method. For each cholesterol concentration, the RDFs calcu- 



lated from the MD simulations were given as an input to the 
IMC procedure to obtain the effective interactions. 

Due to finite size effects, in some cases the RDFs calculated 
from the MD simulations deviate from unity at the cutoff. This 
results in discontinuities in the effective potentials at the cut- 
off. To handle these, we applied a simple shifting scheme 
from 2.0 nm to the cutoff distance to adjust the potentials such 
that they approach zero smoothly at the cutoff. The approach 
used is essentially similar to that presented in Ref. |6J. The 
main difference is that in the present case shifting is applied 
to the potential rather than to the force. 

The IMC does not give reasonable estimates for the effec- 
tive potentials at short interparticle distances where the RDFs 
vanish. In these regions, we replaced the potential given by 
the IMC method by polynomials such that the potential and 
its first derivative are continuous at the edge of the region. Fi- 
nally, the effective potentials were smoothed using the same 
spline-fitting procedure^' as was used for the RDFs to reduce 
statistical noise. We have verified that the potentials are not 
sensitive to the details of the process of obtaining them. Thus 
the above changes can be made without seriously altering the 
resulting RDFs. 

Figure|5] shows the computed effective interaction poten- 
tials. Due to the high level of coarse-graining, the poten- 
tials are soft. The DPPC-DPPC and DPPC-cholesterol poten- 
tials become systematically more repulsive with an increasing 
cholesterol concentration, and the very small attractive com- 
ponent present at low concentrations is lost at higher concen- 
trations. For the cholesterol-cholesterol potentials, the behav- 
ior is more complex. For 29.7 % and 50.0 % concentrations, 
the interaction is mostly repulsive, but for 12.5 % and 20.3 % 
there is weak attraction for r > 0.9 nm up to the cutoff. For 
4.7 %, the interaction is again repulsive for r > 1.1 nm, but 
now there is a weak attraction for small separations. In addi- 
tion, the cholesterol-cholesterol potentials have multiple min- 
ima, whereas the other potentials are much simpler. 



C. Validation 

As any model, the coarse-grained model should be vali- 
dated. This can be done by comparing the results it generates 
to the results from the MD simulations. By construction, the 
CG model should reproduce the short-range structural proper- 
ties of the atomic-scale model. 

Figure |4] shows a comparison between RDFs calculated 
from the MD simulations and those obtained from the CG 
model using canonical MC. The CG simulations contained the 
same number of particles as the MD simulations. We show 
only a few, selected cholesterol concentrations, but for other 
concentrations the results are similar: the agreement between 
the results from MD and CG is excellent at all concentrations, 
as it should. The minor differences near the cutoff arise from 
the use of a shifting function for the potentials. Without the 
shift, the lines coincide up to the cutoff, but in some cases 
there is a small discontinuity in the RDFs at the cutoff. 

We have also calculated the static structure factors com- 
puted over all pairs of particles on a two-dimensional grid. 
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FIG. 4: Comparison between radial distribution functions calculated 
from MD simulations (solid black lines) and from CG model (dashed 
grey lines). Two different cholesterol concentrations are shown: (a) 
12.5 % and (b) 29.7 %. The results for the other concentrations are 
similar. 



To this end, consider a set of particles i — I, . . . ,N whose 
positions are ri . Then the static structure factor defined as 



N N 



^(^) = ^(EEe'^PHfc- {r,-n)]) 



(1) 



.1 j=l 



is given in terms of the reciprocal vector k. The S{k) have 
been calculated for a system whose linear size varies between 
120- 160 nm. This is 24 times the size of the original system 
studied by MD. In all cases the structure factors were found to 
be radially symmetric. The circularly averaged structure fac- 
tors, S{k), for different cholesterol concentrations are shown 
in Fig.|6l These curves are discussed in detail in Sect.|V] At 
this point it is sufficient to note that these calculations confirm 
that the system is isotropic at all concentrations and that there 
is no long-range solid-like order The system is in a fluid-like 
state as it should. 



V. BEHAVIOR AT LARGE LENGTH SCALES 

When the system size is increased from that of the MD 
simulations, several new phenomena are observed. The sim- 
ulations described below were mostly performed on systems 
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FIG. 5: Effective potentials for different pairs of coarse-grained par- 
ticles: (a) DPPC-DPPC, (b) DPPC-cholesterol, and (c) cholesterol- 
cholesterol. 



containing 36 864 particles, corresponding to linear sizes 24 
times those of the original MD simulations. Hence, the linear 
sizes of the systems were 120-160 nm, depending on the con- 
centration. A typical simulation required 50- 100 CPU hours 
on a desktop computer. 

When the system size is increased, there are some minor 
changes in the radial distribution functions. These are illus- 
trated in Fig.0for the 20.3 % cholesterol concentration. For 
other concentrations the results are similar These changes are 
most probably a finite-size effect caused by the small sizes of 
the original atomic-scale systems. It is very likely that if the 
MD simulations could be performed on larger systems, similar 
changes in the RDFs should take place. The figure also shows 
that the RDFs rapidly approach unity at large distances. 

To study possible large-scale organization, we have looked 
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FIG. 6: Total circularly averaged static structure factors computed 
from CG simulations at different cholesterol concentrations. The 
curves are labeled as in Fig.|2| 
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FIG. 7: Changes in radial distribution functions when the size of the 
simulation box is increased at 20.3 % cholesterol. Black solid lines 
correspond to a small system with 64 particles and dashed grey lines 
to a system with a linear size 24 times larger. 



at the static structure factors calculated at different choles- 
terol concentrations. They are shown in Fig. |8] together with 
snapshots of the system at each concentration. The structure 
factors have been calculated over all pairs of molecules, and 
also separately for DPPC-DPPC, cholesterol-cholesterol, and 
DPPC-cholesterol pairs. To ensure that the static structure fac- 
tors do not depend on the initial configuration, we have rerun 
the calculations using several different initial states. We find 
that different initial configurations lead to results that are con- 
sistent with each other 

At 12.5 % and 20.3 % cholesterol, the snapshots in Fig. |8l 
suggest that there are domains where the local concentration 
of cholesterol is higher than in other regions. At these con- 
centrations the cholesterol-cholesterol structure factor shows 
a rather wide peak at small k, supporting our interpretation 
of the presence of cholesterol-rich and cholesterol-poor do- 



mains. The maximum of the peak corresponds to length scales 
of the order of 20 nm or more. A more precise analysis is dif- 
ficult, since even for the largest systems we have studied, with 
linear sizes of approximately 280 nm, the peak is rather broad, 
and the small k side of the peak is not fully clear due to fluc- 
tuations at large length scales. 

We have, however, studied how the peak behaves if the sys- 
tem size is varied by calculating the static structure factors for 
four systems with linear sizes 6, 12, 24 and 48 times that of 
the original MD simulations. In this finite size scaling analy- 
sis, we found no clear change in either the position or shape 
of the peak. 

The formation of cholesterol-rich and cholesterol-poor do- 
mains at intermediate cholesterol concentrations is in agree- 
ment with the phase diagram, see Fig. ^ According to the 
phase diagram we should, at T = 323 K, expect both the MD 
and the CG model to display coexisting Id and lo phases. In 
the present case, we cannot directly distinguish between the 
Id and lo phases, since we have not included the ordering of 
the DPPC tails in the model. However, we might argue as fol- 
lows to establish that the cholesterol-rich phases should be lo, 
while the cholesterol-poor must be Id. The ordering effect of 
cholesterol on the phospholipid tails has been clearly demon- 
strated: the higher the cholesterol concentration, the more or- 
dered the tails. '^-^^ It is thus plausible that the tails should be 
more ordered in the cholesterol-rich regions. 

To study domain formation in more detail, we have com- 
puted probability distribution functions for finding square do- 
mains with a linear size t. Several different system sizes t 
were considered. If there were domains of clearly differ- 
ent compositions, the distribution should display two peaks. 
These studies do not provide direct support for formation of 
domains (data not shown). At each concentration, and for 
each i, the distribution is very close to a Gaussian. However, 
with 12.5 % and 20.3 % cholesterol the variance of the distri- 
bution is clearly larger than for 4.7 % or 29.7 % cholesterol. 
Thus, it is possible that the observed domains are caused by 
large-scale fluctuations in the cholesterol concentration. Fur- 
thermore, at 12.5 % the distribution is not quite symmetric for 
I < 20 nm. This may be due to two peaks that are located 
very close to each other. 

We have also investigated the local concentrations of dif- 
ferent types of molecules in the vicinity of a given type of 
molecule as in a study by de Vries et al.^^ More specifically, 
we have studied the respective numbers of different types 
of molecules among the n nearest neighbors of DPPC and 
cholesterol molecules, for different values of n. The average 
number of molecules of type A among the n nearest neigh- 
bors of a molecule of type B was compared to the value in a 
purely random configuration. In all cases the average number 
of cholesterol molecules in the vicinity of another cholesterol 
molecule was somewhat smaller than it should be in a random 
configuration. For n = 6, there were no significant differ- 
ences between different cholesterol concentrations. For larger 
n, say n = 15 or n = 30, the difference between the value 
from an actual simulation and the value from a random con- 
figuration was clearly smaller at the concentrations where we 
observe cholesterol-rich and cholesterol-poor domains than at 



^ 
^ 



i.Q 
O.S 
Q.6 
0.4 

M 

Q.Q 



' I ' I ' I ■ I ' I ' I ' I 



■ 1 .■Sib'^--1-- 1 "'I 
«' ^ . i/s, Ti. , •■. -I 




A' [nm ] 



it [nm' ] 



FIG. 8: Static structure factors calculated for different cholesterol concentrations: (a) 4.7 %, (b) 12.5 %, (c) 20.3 %, and (d) 29.7 %. In addition 
to the total structure factor computed over all pairs of particles (solid line), structure factors calculated over DPPC-DPPC (dotted), cholesterol- 
cholesterol (dashed) and DPPC-cholesterol (dash-dotted) pairs are shown. Also a snapshot of the system is shown for each concentration. In 
the snapshots, each cholesterol molecule is represented as a single dot, while DPPC molecules are not shown. 



Other concentrations. This observation may be interpreted as 
support to the existence of domains. The reason why the dif- 
ference is not visible in the case of n = 6 is probably that such 
a small neighborhood represents the area from which other 
cholesterol molecules are largely excluded by the repulsive 
interactions between the cholesterol molecules. 

Based mainly on fluorescence measurements, it has been 
proposed that at certain cholesterol concentrations the choles- 
terol molecules could adopt a regular arrangement in parts 
of the bilayer''^ We have studied whether or not our coarse- 
grained model shows such organization for a few of the 
proposed "magical" concentrations. The nearest-neighbor 
analysis described above suggests that finding a choles- 
terol molecule in the immediate vicinity of another choles- 
terol molecule is lower compared to a random configuration, 
whereas the probability of finding a cholesterol molecule next 
to a DPPC molecule is higher. Such a situation would oc- 
cur if there were superlattices in the system. Nevertheless, 
we have seen no evidence for any regular lattice-like ordering 
of cholesterols. Despite this conclusion, we are looking for- 
ward to further studies using other models. In our model the 



DPPC particles are radially symmetric, while in reality the 
two hydrocarbon tails of the molecules may be important for 
the occurrence of superlattices, if they do exist. 



In the case of 50 % cholesterol the situation is more com- 
plex. Visual inspection of snapshots clearly indicates that 
the DPPC and cholesterol molecules phase separate. All 
other studies we have performed support this view. There 
are at least two alternative explanations for this phenomenon. 
First, it is possible that the atomistic MD simulations do 
not, at this high cholesterol concentration characterized by a 
very small lateral diffusion coefficient, adequately sample the 
phase space. If this is the case, it may lead to errors in the 
RDFs extracted from the atomic-level simulations. As a con- 
sequence, the effective potentials given by the IMC method 
would not describe the true behavior On the other hand, 
there is experimental evidence for the formation of crystalline 
cholesterol domains at very high cholesterol concentrations.— 
Also this could be the reason for the observed phase separa- 
tion. Based on the current simulations, it is difficult to say 
which of these, if indeed any, is the case. 
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VI. DISCUSSION 

The coarse-graining approach described here has several 
advantages. It is a systematic method to generate a mesoscale 
model that is linked to atomic-scale information. It is ad- 
justable, as it allows the user to control the level of coarse- 
graining and the number of degrees of freedom to be included 
in the model. The approach is general in the sense that it can 
be applied to a wide range of systems. The only prerequisite is 
that we should somehow acquire radial distribution functions 
for pairs of the coarse-grained particles. These are available 
from atomic-scale simulations, but in certain cases also exper- 
iments could be used to derive the RDFs. 

The most important advantage, however, is speed. Monte 
Carlo simulations using the coarse-grained model are several 
orders of magnitude faster than MD simulations of the orig- 
inal atomic-scale model. The speedup can be estimated by 
considering the decay time of the fluctuations in the area of a 
single molecule, where the area is calculated, through Voronoi 
analysis. ^"'^^ For the atomic-scale model, this decay time is 
approximately 0.7 ns,^^ whereas in the case of the CG model 
similar decay is observed after one MC step. The CPU time 
required for generating a 0.7 ns trajectory for the atomic-scale 
system using MD is around 70 h, whereas for a CG model of 
the same linear size as the atomic-scale system, approximately 
25 000 Monte Carlo steps can be taken during a minute. Thus 
the speedup can be estimated to be eight orders of magnitude. 

Despite all advantages, there are limitations. An impor- 
tant point is that any problems in the atomic-scale molecu- 
lar dynamics simulations are transferred to the coarse-grained 
model. In particular, if there is any artificial ordering in the 
atomic-scale simulations, either due to a small system size, 
poor sampling of the phase space, or sloppy treatment of elec- 
trostatic interactions,^ ''^^ the radial distribution functions will 
be erroneous. As a consequence, the effective interaction po- 
tentials will be affected. 

Another limitation may be the concentration-dependence of 
the effective potentials, see Fig. |5] It is not trivial in what 
concentration ranges the effective potentials are valid. The 
worst case scenario is that when the cholesterol concentration 
is altered ever so slightly, the interaction potentials must be 
rederived starting from time-consuming atomic-scale molec- 
ular dynamics simulations. We have investigated whether the 
effective potentials can be used at concentrations other than 
those at which they were determined. The results suggest 
that they may be used at nearby concentrations, but if phase 
boundaries are crossed, problems will arise. For instance, 
when using the effective potentials determined for the system 
with 29.7 % cholesterol, a system in the lo phase, to simu- 
late a system at 20.3 % cholesterol, which should be in the 
coexistence region, no long-range structure appears. When, 
on the other hand, potentials determined at 20.3 % cholesterol 
are used with 29.7 % cholesterol, the long-range structure is 
still present at the higher concentration, although the peak is 
more shallow. The situation is similar when we compare the 
4.7% and 12.5% concentrations. When the potentials deter- 
mined at 12.5 % are used for 20.3 % cholesterol or vice versa, 
i. e., both concentrations are in the coexistence region, we still 



observe domains. The detailed form of the static structure fac- 
tor will, however, be altered. We may therefore conclude that 
the effective potentials cannot be used for mapping the precise 
phase boundaries of a given system. Similar conclusions can 
be drawn from a study of the temperature dependence of the 
effective interaction potentials. 

An additional problem is the implementation of dynam- 
ics. In this study, we considered only structural quantities, 
and needed not to worry about the choice of dynamics. To 
study lateral diffusion it is necessary to incorporate realistic 
dynamics into the system. It is well-known that this can be 
very difficult in the case of coarse-grained models. In addi- 
tion to the MC studies reported here, we attempted to study 
lateral diffusion and include dynamics into our model. In 
choosing the dynamics, the following criteria should be met. 
First, as explained in Sect. |III| the dissipative nature of the 
system should be taken into account. In addition, the simula- 
tions should be run in the canonical ensemble. We tested sev- 
eral thermostats that comply to the above requirements: stan- 
dard and generalized Brownian dynamics, ^^ and an approach 
similar to Andersen's thermostat for temperature coupling.— 
None of these could be tuned to give realistic dynamics for 
the system with rattling-in-the-cage movement and separate 
jump events. Both the standard Brownian dynamics approach 
and the Andersen scheme require very large friction or cou- 
pling parameters to give diffusion coefficients of the same or- 
der of magnitude as those obtained from the molecular dy- 
namics simulations, or alternatively, to match the decay times 
of velocity autocorrelation functions. Such high values of the 
parameters completely determine the dynamics at short time 
scales, and the interactions between the particles only give 
rise to small corrections at long time scales. The general- 
ized Brownian dynamics can be used to force the short-time 
dynamics to match that of molecular dynamics simulations, 
but this does not mend the problem of unrealistic dynamics 
at longer time scales. We thus found no simple way of im- 
plementing realistic dynamics into the coarse-grained model. 
Reducing the degree of coarse-graining and introducing more 
detail into the model, e.g., by including the tails of DPPC 
and cholesterol molecules, might help to solve the problem. 
The presence of tails would greatly increase the friction be- 
tween the molecules, and possibly allow for entanglements. 
Both effects would help in slowing down the unrealistically 
fast dynamics. 

The results motivate further studies of the coarse-graining 
approach. There are several possible directions in which the 
model could be developed. Several of these could be real- 
ized without additional MD simulations. A possible line of 
development is the inclusion of the conformational degrees of 
freedom of the DPPC molecules in the model. This could be 
done in the spirit of the model by Nielsen et al.,^"^ i. e., by giv- 
ing each DPPC molecule two possible states: an ordered and 
a disordered state. We would have three kinds of particles and 
a total of six pairwise potentials to determine. Another possi- 
ble modification would be to include the two tails of the lipid 
molecules as separate particles and possibly model the head 
group as a third particle. Results from simulations of such 
models could be compared to the present study to assess the 
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possible benefits of such modifications, and to gain further in- 
sight into the coarse-graining process. Additionally, to better 
understand the underlying reasons for domain formation, it is 
natural to ask what the relative roles of entropic and energetic 
contributions in this process are. Since the close-packed areas 
of the CG particles are not well defined, a study of this broad 
and non-trivial issue is beyond the scope of this work and will 
be discussed elsewhere. 

The coarse-graining approach presented here could also be 
applied to other lipid bilayer systems to compare their be- 
havior to the DPPC/ cholesterol bilayer A particularly inter- 
esting system is the sphingomyelin (SM) / cholesterol bilayer, 
and ultimately the PC /SM/ cholesterol ternary mixture. The 
study of these systems at large length scales would be par- 
ticularly interesting, because experimental results suggest the 
existence of specific interactions between SM and cholesterol 
moleculesi2£ These interactions are believed to forward the 
formation of domains and lipid rafts at high concentrations of 
SM and cholesterol. This line of development is limited by the 
computational cost of the underlying MD simulations. Cur- 
rently this confines us to study quite simple bilayer systems. 
With increasing computer power and improvements in exist- 
ing simulation methods, studies of ternary mixtures of lipids, 
as well as studies of systems containing membrane proteins, 
are becoming feasible. 



VII. SUMMARY AND CONCLUSIONS 

We have used a systematic approach for constructing 
coarse-grained models for DPPC / cholesterol bilayers. The 
central ingredient is the application of the Inverse Monte Carlo 
method, which can be used for finding effective interactions 
such that the coarse-grained model reproduces given radial 
distribution functions. The approach allows easy tuning of 



the level of coarse-graining, and it can be applied to a wide 
range of systems. 

The radial distribution functions given as input to the IMC 
method have been extracted from detailed atomic-level molec- 
ular dynamics simulations. The effective interactions found 
using the IMC method can then be used to simulate the system 
on much longer length scales. We have found that the coarse- 
grained model thus constructed is in favor of the formation of 
cholesterol-rich and cholesterol-poor domains at intermediate 
cholesterol concentrations, in agreement with the phase dia- 
gram of the system. We have also explored the limitations of 
the constructed coarse-grained model. 

As for further studies, it would be interesting to see how 
modifications such as the inclusion of the conformational de- 
grees of freedom of the DPPC tails would influence the behav- 
ior of the model. Similar models could also be constructed for 
other many-component lipid bilayer systems, and their behav- 
ior compared to the present study. This could yield valuable 
new information on both the systems under study and the suit- 
ability of the IMC method for coarse-graining in general. 
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